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ABSTRACT 

In  certain  geophysical  contexts  such  as  lava  lakes  and  mantle  convection,  a  cold, 
viscous  boundary  layer  forms  over  a  deep  pool.  The  following  model  problem  investigates 
the  buoyant  instability  of  the  layer.  Beneath  a  shear-free  horizontal  boundary,  a  thin  layer 
(thickness  d\)  of  very  viscous  fluid  overlies  a  deep  layer  of  less  dense,  much  less  viscous 
fluid;  inertia  and  surface  tension  are  negligible.  After  the  initial  unstable  equilibrium 
is  perturbed,  a  long-wave  analysis  describes  the  growth  of  the  disturbance,  including  the 
nonlinear  effects  of  large  amplitude.  The  results  show  that  nonlinear  effects  greatly  enhance 
growth,  so  that  initial  local  maxima  in  the  thickness  of  the  viscous  film  grow  to  infinite 
thickness  in  finite  time,  with  a  timescale  8/i/ Ap gd\.  In  the  final  catastrophic  growth  the 
peak  thickness  is  inversely  proportional  to  the  remaining  time.  (A  parallel  analysis  for 
fluids  with  power-law  rheology  shows  similar  catastrophic  growth.)  While  the  small-slope 
approximation  must  fail  before  this  singular  time,  the  failure  is  only  local,  and  a  similarity 
solution  describes  how  the  peaks  become  downwelling  plumes  as  the  viscous  film  drains 
away. 
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BUOYANT  INSTABILITY  OF  A  VISCOUS  FILM  OVER  A  PASSIVE  FLUID 
1.  INTRODUCTION 

This  work  examines  the  strongly  nonlinear  effects  of  finite  amplitude  in  the  Rayleigh- 
Taylor  instability  of  a  horizontal  viscous  film  under  a  shear-free  boundary  and  over  a  much 
less  viscous  fluid.  Inertia  and  surface  tension  are  neglected,  and,  in  the  parameter  range 
considered,  the  motion  is  limited  by  normal  stresses  in  the  more  viscous  fluid.  The  analysis 
exploits  the  fact  that  the  most  unstable  wavelengths  are  long  compared  to  the  thickness 
of  the  film.  The  results  show  how  the  growth  of  disturbances  to  the  interface  becomes 
greatly  enhanced  when  the  disturbance  amplitude  becomes  large,  leading  to  the  formation 
of  downwelling  sheets  or  plumes  in  a  finite  time. 

The  motivation  for  this  problem  comes  from  certain  geophysical  situations,  partic¬ 
ularly  the  stability  of  the  Earth’s  lithospheric  plates.  In  simplified  terms,  the  oceanic 
lithosphere  (tectonic  plates)  can  be  considered  a  cold,  stiff  thermal  boundary  layer  above 
the  convecting  mantle.  Where  two  plates  come  together,  one  subducts  under  the  other 
and  flows  downward  due  to  its  negative  buoyancy.  The  question  of  how  a  new  subduction 
zone  is  formed,  how  one  large  plate  may  break  into  two  and  thus  allow  some  of  the  dense 
material  to  flow  back  down  into  the  mantle,  is  not  yet  resolved.  Other  closely  related 
geophysical  situations  include  the  surfaces  of  lava  lakes,  thermal  convection  in  the  mantles 
of  other  planets,  and  possibly  convection  in  the  Earth’s  solid  core. 

This  work  examines  a  simple  model  of  one  possible  mechanism  for  the  initiation  of 
subduction:  the  Rayleigh-Taylor  instability.  In  this  model,  the  lithosphere  and  the  mantle 
are  treated  as  distinct,  highly  viscous  fluids,  the  lithosphere  being  denser  (and  much  more 
viscous)  than  the  mantle.  In  this  unstable  configuration,  any  variations  in  the  lithosphere 
thickness  tend  to  grow,  and  when  the  thickness  variations  become  significant,  nonlinear  ef¬ 
fects  cause  the  growth  to  accelerate  catastrophically,  giving  downwelling  regions  (modelling 
subduction). 

A  more  realistic  model  of  the  lithosphere  may  be  a  fluid  with  a  power-law  rheology, 
with  most  of  the  layer  being  very  stiff,  but  yielding  more  readily  in  regions  of  rapid  de¬ 
formation.  A  long-wave  analysis  for  this  case  (Appendix  B)  shows  that  the  growth  in  the 
deforming  regions  again  becomes  catastrophic  due  to  finite-amplitude  effects,  as  in  the 
Newtonian  case. 

A  companion  work  (Canright  and  Morris,  in  prep.)  considers  a  different  model  for 
the  lithosphere:  as  a  thermal  boundary  layer  growing  under  a  suddenly  cooled  horizontal 
boundary,  where  the  fluid  viscosity  depends  strongly  on  temperature.  The  long-wave  anal¬ 
ysis  shows  that,  again,  the  nonlinear  effects  of  finite  amplitude  give  catastrophic  growth, 
yielding  sheets  in  finite  time.  For  that  case  the  force  balance  and  the  resulting  growth  of 
peaks  is  essentially  the  same  as  that  considered  here,  because  thermal  diffusion  becomes 
unimportant  where  the  layer  is  thick.  The  Rayleigh-Taylor  problem  considered  here  is  the 
simplest  example  of  this  dynamic  balance. 

The  instability  of  a  dense  fluid  supported  by  a  lighter  one  has  been  studied  extensively, 
beginning  with  the  analysis  of  Rayleigh  (1883).  Taylor  (1950)  noted  that  acceleration 
would  give  the  same  effect  as  gravity,  in  agreement  with  the  experiments  of  Lewis  (1950). 
A  thorough  introduction  to  the  linear  theory  is  given  by  Chandrasekhar  (1961,  Ch.  X), 
and  recent  surveys  of  previous  work  are  given  by  Sharp  (1984)  and  Kull  (1991).  Many 
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of  the  analyses  (e.g.,  Bellman  and  Pennington,  1954,  and  Menikoff  et  al.,  1978)  have 
focused  on  the  initial  small-amplitude  growth,  where  linearized  equations  are  appropriate. 
One  approach  for  examining  the  nonlinear  effects  of  large  amplitude  is  direct  numerical 
simulation,  e.g.,  Harlow  and  Welch  (1966)  and  Tryggvason  (1988),  or  through  boundary- 
integral  computations,  e.g.,  Baker  and  Meiron  (1984),  and  Newhouse  and  Pozrikidis  (1990). 

To  investigate  the  nonlinear  effects  analytically,  the  most  common  approach  has  been 
the  perturbation  method,  where  various  quantities  are  expressed  as  power  series  in  the 
small  initial  amplitude  or  slope,  giving  linearized  equations  for  each  order  (e.g.,  Emmons 
et  ah,  1960).  However,  a  wide  variety  of  other  methods  have  been  employed  (primarily 
for  inviscid  fluids)  e.g.,  least-squares  approximation  (Kull,  1986),  averaging  (Drazin,  1969), 
multiple  time  scales  (Nayfeh,  1969),  strained  coordinates  (Amaranath  and  Rajappa,  1976), 
generalized  coordinates  (Dienes,  1978),  Lagrangian  formulations  (Ott,  1972),  and  heuristic 
models  (Baker  and  Freeman,  1981,  Aref  and  Tryggvason,  1989).  Also,  Dussan  V.  (1975) 
used  energy  methods  to  determine  the  stability  of  disturbances  of  arbitrary  amplitude. 
In  the  above  studies,  inertia  is  important.  The  situation  we  consider,  of  a  thin  viscous 
layer  in  creeping  flow,  has  been  examined  primarily  in  the  geophysics  literature,  e.g.,  the 
experiments  and  weakly  nonlinear  theory  of  Whitehead  and  Luther  (1975). 

Our  analysis  is  equivalent  to  a  perturbation  expansion  in  the  interfacial  slope,  which 
for  long  waves  remains  small  even  for  large  amplitudes.  But  unlike  the  standard  approach, 
the  interface  position  is  not  expressed  as  a  series;  hence  the  (leading-order)  interface  con¬ 
ditions  are  nonlinear,  being  applied  at  the  moving  interface  rather  than  expanded  about 
a  plane.  In  this  way  the  nonlinear  effects  of  finite  amplitude  are  explicitly  included.  We 
derive  the  small-slope  equations  for  three  dimensions,  and  give  solutions  for  the  two- 
dimensional  and  axisymmetric  cases.  The  results  show  that  thickness  maxima  grow  to 
infinite  thickness  in  finite  time;  thereafter  the  layer  drains  through  these  downwelling  re¬ 
gions.  (Of  course,  the  small-slope  assumption  must  fail  before  this,  but  the  failure  is 
localized  in  asymptotically  narrow  regions  with  negligible  effects  on  the  dynamics  in  the 
rest  of  the  layer.) 

2.  PROBLEM  STATEMENT 

Two  horizontal  layers  of  distinct  Newtonian  incompressible  fluids  are  bounded  above 
and  below  by  horizontal  shear-free  boundaries  (see  figure  1).  The  upper  fluid  (of  density 
Pi,  viscosity  p,\.  and  average  layer  thickness  cli)  is  denser  and  much  more  viscous  than 
the  lower  fluid  (of  density  p2  <  p\  and  viscosity  p2  -C  /  / 1 ) ,  and  the  upper  layer  is  very 
thin  (c/i  -C  d2).  Surface  tension  is  neglected,  and  both  fluids  are  assumed  to  be  so  viscous 
that  we  can  neglect  inertia.  The  former  assumption  requires  that  yfc2 / Apg  <C  1,  while  the 
latter  requires  ApgL3  /  p2u2  <A  1,  where  7  is  the  surface  tension,  k  is  the  wavenumber,  L  is 
the  largest  length  scale,  A p  =  pi  —  p2,  g  is  the  acceleration  of  gravity,  and  v  is  the  smaller 
kinematic  viscosity.  (The  latter  restriction  is  easily  satisfied  in  mantle  flow.) 

Initially  both  fluids  are  at  rest  in  unstable  equilibrium.  At  t  =  0  the  interface  is  slightly 
disturbed,  after  which  the  position  of  the  interface  is  d(x,y,t).  (Long-wave  solutions  will 
be  given  only  for  the  special  cases  of  two-dimensional  and  axisymmetric  disturbances.) 

A  reduced  pressure  p  is  defined  in  terms  of  the  total  pressure  P  by: 

P(x,y,z,t)=p(x,y,z,t)  +  p2gz  (2.1) 
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so  that  in  equilibrium  p  is  a  constant  in  fluid  2.  Then  the  governing  equations  are: 


in  fluid  1:  V(p  +  Apgz)  =  /ii  V2u  (2.2a) 

in  fluid  2:  Vp  =  /i2V2u  (2.2b) 

V  u  =  0  (2.3) 

where  u  is  the  velocity  vector  with  components  (u.  v.  w).  The  boundaries  exert  no  shear, 
and  across  the  interface,  which  moves  as  a  material  surface,  both  velocity  and  stress  are 
continuous,  hence: 


at  z  =  0  and  at  z  =  d\  +  d,2  :  w  —  uz  —  vz  =  0  (2.4a) 

at  z  =  S(x,  y,  t)  :  [u]  |  ^  =  0 ,  (2.4b) 

[cr.j-nj]  |^  =  0,  (2.4c) 

5t  +  u5x  +  vSy  =  w  (2.4d) 


2 

where  subscripts  indicate  partial  derivatives,  the  brackets  []  |~  indicate  the  change  in  the 
enclosed  quantity  across  the  interface  from  fluid  1  to  fluid  2,  at]  is  the  reduced  stress 
tensor  (using  the  reduced  pressure  p),  and  n  is  a  unit  vector  normal  to  the  interface,  with 
components  n% .  (If  the  upper  boundary  is  instead  a  free  surface  with  no  shear  and  zero 
pressure,  then  in  the  long-wave  analysis  below,  the  surface  deflection  is  proportional  to  the 
two-fluid  interface  deflection  in  the  ratio  A  p/p,  which  is  assumed  to  be  small.  Also  note 
that  the  boundary  condition  at  z  —  d\  +  (h  is  irrelevant  to  the  long-wave  analysis  below, 
provided  fluid  2  is  sufficiently  deep.)  Specification  of  the  initial  interface  position  S(x,  y.  0) 
completes  the  problem  definition. 

3.  LONG- WAVE  ANALYSIS 

The  disturbance  can  grow  in  a  variety  of  different  ways,  depending  on  the  viscosity 
ratio,  the  depth  ratio,  and  the  dimensionless  wavenumber: 

a  =  p  2/ Pi  ,  /3  =  c/2/di  ,  e  =  kdi  (3-1) 

where  the  disturbance  has  a  characteristic  wavenumber  k.  Here  we  focus  on  the  case  where 
the  lower  fluid  layer  is  much  less  viscous  and  deeper  than  the  upper  (a  <C  1,  j3  V>  1). 

The  linearized  small-amplitude  solution  given  in  Appendix  A  shows  that  for  this  case 
the  fastest  growing  wavelength  is  long  compared  to  d\ .  In  fact,  the  linearized  growth  rate 
is  nearly  constant  over  a  broad  range  of  wavelengths: 

max  (a,  y/a/(3)  <  e  «  1  (3.2) 

as  illustrated  in  figure  7a.  In  this  range,  the  growth  is  limited  by  normal  stresses  in  the 
upper  fluid,  which  moves  nearly  horizontally,  while  the  lower  fluid  is  passively  moved  by 
the  interface.  Outside  this  range,  for  waves  short  compared  to  di,  the  growth  is  reduced 
because  only  a  fraction  of  layer  1  is  mobilized,  and  for  long  enough  waves  the  viscous 
resistance  of  fluid  2  slows  the  growth. 
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We  examine  the  finite-amplitude  growth  of  long-wave  disturbances  in  this  fastest- 
growing  range,  exploiting  the  fact  the  slope  of  the  interface  remains  small  even  for  large 
amplitudes  (until  the  disturbance  grows  to  the  order  of  the  initial  wavelength).  What 
follows  is  the  leading-order  asymptotic  analysis,  containing  only  the  dominant  balance  of 
forces.  This  balance  is  insensitive  to  wavelength  (the  linearized  growth  rate  is  constant- 
to  within  0(e4)),  so  this  analysis  cannot  predict  the  single  wavelength  giving  maximum 
growth. 

The  small  slope  of  the  interface  greatly  simplifies  the  dynamics  of  the  thin  upper  layer, 
which  controls  the  growth.  In  fluid  1,  the  horizontal  length  scale  (/c_1)  is  much  larger  than 
the  vertical  (di),  so  by  continuity  the  horizontal  velocity  components  (scale  U ,  say)  are 
large  relative  to  the  vertical  components  (eU).  Then  the  momentum  equation  implies  the 
vertical  variation  of  reduced  pressure  is  negligible  (0(e))  compared  to  horizontal  variations. 

In  fluid  2,  the  horizontal  velocity  and  length  scales  are  those  of  the  interface  ( U  and 
/c-1);  the  vertical  length  scale  is  comparable  (except  for  waves  long  compared  to  da-  where 
the  depth  is  the  vertical  scale).  Then  the  shear  stress  from  fluid  2  scales  as  p2 uz  ~  fi2 kU 
(or  ~  ^20/^2  in  the  latter  case).  Since  the  shear  in  fluid  1  scales  as  pi uz  ~  piwx  ~  pi ekU, 
then  the  shear  from  fluid  2  will  be  negligible  provided  a  • C  e  (or  a/ (3  -C  e2  if  (3e  <  1), 
which  is  just  the  lower  bound  on  wavenumber  assumed  above.  For  wavelengths  in  this 
range,  the  interface  motion  is  controlled  by  the  dynamics  of  the  upper  fluid,  for  which  the 
lower  fluid  is  effectively  passive  and  hydrostatic. 

Then  in  layer  1,  both  the  upper  boundary  and  the  interface  appear  shear-free,  so  the 
horizontal  velocity  components  are  independent  of  z  (to  0(e2)).  as  are  the  components 
rxx.  rxy,  Tyy,  rzz  of  t he  deviatoric  stress  tensor  t-13  .  (This  scaling  analysis  is  extended  to 
other  wavelengths  in  Canright,  1987.) 

While  the  long-wave  equations  can  be  derived  by  expanding  the  velocity  and  stress  in 
powers  of  e,  the  meaning  of  the  result  is  better  understood  as  a  force  balance  on  a  small 
column  of  layer  1,  as  shown  in  figure  2.  Here  only  horizontal  forces  are  considered,  and 
as  before  the  total  stress  everywhere  has  been  reduced  by  the  hydrostatic  pressure  of  fluid 
2,  which  doesn’t  effect  the  horizontal  balance.  Then  there  are  no  forces  acting  on  the  top 
and  bottom  of  the  column,  and  so  the  net  force  acting  on  the  sides  must  be  zero.  The 
differential  form  of  this  balance  can  be  written  in  terms  of  a  two-dimensional  reduced  force 
tensor  F,3  which  results  from  integrating  the  reduced  stress  tensor  across  the  layer 

Fij=  [  [~(p  + Apgz)Sij +  Tij]  dz  ,  i,j  =  1,2  (3.3) 

Jo 

where  Sl3  is  the  Kronecker  delta.  The  normal  stress  condition  at  the  interface 


at  z  —  5  :  —  (p  +  Apg5)  +  tzz  =  0  (3-4) 

can  be  used  to  eliminate  p,  so  the  reduced  force  is 

Fi3  =  5(-P5i:j  +  Tij )  (3.5a) 

P  =  -Apg5/2  +  tzz  (3.5b) 
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where  P  is  the  average  pressure  and  rzz  is  evaluated  at  the  interface.  Then  the  horizontal 
force  balance  means  the  reduced  force  tensor  has  zero  divergence 


Fij,j  =  0  ,  i,j  =  1,2 

The  kinematic  interface  condition  becomes,  on  eliminating  w  using  continuity: 


(3.6) 


5t  +  (du)x  +  (Sv)y  —  0 


(3.7) 


These  two  equations  (along  with  conditions  at  the  edges  of  the  layer)  govern  the  growth  of 
long-wave  disturbances,  including  the  nonlinear  effects  of  finite  amplitude,  with  a  relative 
error  of  0(e).  This  conclusion  is  independent  of  rheology,  and  relies  only  on  the  small- 
slope  approximation  (valid  until  amplitudes  are  comparable  to  the  wavelength)  and  the 
assumption  that  fluid  2  is  dynamically  unimportant. 

The  dimensionless  form  of  (3.6)  becomes,  for  Newtonian  rheology, 

[S2  +  5{ux  +  \vy)\  x  +  \  [S(uy  +  vx)]y  =  0  (3.8a) 

[$2  +  5(vy  +  \ux)\  y  +  \  [5(uy  +  vx)]x  =  0  (3.8b) 

while  (3.7)  remains  the  same  under  non-dimensionalization.  (A  fluid  with  power-law  rhe¬ 
ology  is  considered  in  Appendix  B.)  Here  we  used  a  vertical  length  scale  of  di,  a  time  scale 
of  8f.i/ Apgdi  (corresponding  to  half  the  small- amplitude  growth  rate),  and  any  horizontal 
length  scale.  In  this  simplified  problem,  the  layer  thickness  6  is  analogous  to  variable  den¬ 
sity,  as  shown  by  (3.7),  but  also  acts  somewhat  like  variable  viscosity,  as  well  as  a  body 
force,  iu  (3.8). 

If  the  initial  disturbance  to  the  interface  is  axisymmetric  or  varies  only  in  one  direction, 
the  dimensionality  of  the  problem  is  further  reduced,  allowing  simpler  solutions  to  illustrate 
the  nonlinear  dynamics.  The  axisymmetric  equations  are 

(3.9) 
(3.10) 

where  u  is  the  radial  velocity  component  and  r  the  radial  coordinate. 

The  long-wave  versions  of  the  two-dimensional  equations 


52  +  5 


N, 

r 

(■ rSu ) 


=  -5  - 

—  2Ur  j. 


=  o 


(52+5ux)x  =  0  (3.11) 

St  +  (Su)x  =  0  (3.12) 


allow  greater  simplification.  The  first  shows  that  the  reduced  force  (scalar)  F  =  52  +  5ux 
acting  on  any  vertical  cross-section  of  the  layer  is  independent  of  position,  though  it  may 
depend  on  time;  the  value  of  F(t)  depends  on  the  conditions  at  the  ends  of  the  layer. 

We  can  eliminate  u  by  adopting  a  Lagrangian  formulation,  where  the  fluid  “particle” 
in  this  case  is  a  material  cross  section  of  the  layer.  Then  the  new  variables  are  (xo,t), 
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where  Xo  refers  to  a  particular  fluid  section  by  its  initial  position.  The  resulting  ordinary 
differential  equation  describes  the  growth  of  a  moving  cross  section: 

hl=«S2-^(t)  (3.13) 

where  =  St  +  u5x  and  de(t)  =  \J F(t)  is  the  current  equilibrium  thickness,  i.e. ,  that 
with  no  tendency  to  grow  or  shrink.  This  result  shows  that  the  vertical  motion  of  the 
interface  depends  only  on  the  local  layer  thickness  and  the  current  equilibrium  thickness 
(which  may  depend  on  the  overall  shape  of  the  whole  layer).  Thus  any  fluid  cross  section 
does  not  care  what  its  immediate  neighbors  are  doing,  and  the  growth  is  insensitive  to 
wavelength,  as  expected. 

The  horizontal  motion  of  a  cross  section  depends  on  the  state  of  the  whole  layer. 
The  position  of  each  fluid  cross  section  can  be  followed  by  integrating  the  strain.  By 
conservation  of  mass,  5odxo  —  5dx,  where  $0(^0)  =  d(xo,t  —  0)  is  the  initial  thickness  of 
the  section,  so  the  strain  is  dx/dx o  —  5q/5  and  the  position  is  given  by: 


x(x0,  t ) 


M£o) 


dio 


(3.14) 


where  the  x  origin  is  chosen  at  some  stationary  fluid  cross  section. 


4.  SOLUTIONS 

To  show  most  clearly  the  effects  of  finite  amplitude,  we  concentrate  on  the  simplest 
situation,  where  the  initial  disturbance  depends  only  on  x.  We  give  solutions  to  (3.13)  for 
three  types  of  conditions  at  the  ends  of  the  layer,  for  which  de(t)  is  simple  to  evaluate.  For 
the  first,  the  layer  is  of  infinite  extent  but  the  disturbance  is  localized,  so  that  far  away 
the  layer  remains  in  (unstable)  equilibrium,  and  de  —  1.  Next  is  the  case  of  a  periodic 
disturbance;  then  de(t)  is  an  integral  property  of  the  shape  of  the  layer.  The  third  case 
treats  a  layer  of  finite  extent,  with  abrupt  ends  surrounded  by  fluid  2,  so  de  —  0;  this 
solution  also  applies  approximately  whenever  6  grows  large  enough  that  de  is  negligible. 
Lastly,  to  illustrate  one  significant  difference  when  the  disturbance  varies  in  more  than  one 
direction,  we  show  a  solution  to  (3.9,  3.10),  the  axisymmetric  case. 


4.1.  LOCALIZED  DISTURBANCE 


The  most  illustrative  case  is  an  infinite  layer  with  a  disturbance  of  finite  extent.  Then 
the  force  on  the  ends  of  the  disturbed  region  remains  constant:  de  =  1.  This  gives  a  simple 
solution: 


S{x0,t) 


8o(x0)  —  tanh(t) 

1  —  $0(^0)  tanh(t) 


(4.1) 


where  again  is  the  initial  thickness  and  Xq  the  initial  position  of  a  fluid  section.  This 
solution  predicts  a  singularity  in  finite  time:  at  the  critical  time  f *  =  -f  In  |  (<50  +  1)/ (80  —  1)  I 
the  thickness  5  goes  to  00  or  0,  depending  on  whether  the  initial  thickness  40  was  greater  or 
less  than  1.  (Of  course,  the  small-slope  approximation  must  fail  when  the  peak  thickness 
gets  large  enough;  this  point  is  addressed  in  the  next  section.) 
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The  early  growth  of  a  disturbance  of  small  amplitude  a(x o,  t)  =  <5(xo,  t)  —  1  is  roughly: 

a(x0,  t )  «  a0(xo)e2t[l  +  a0(x0)(e2t  -  l)/2]  (4.2) 

Initially  this  gives  the  exponential  growth  of  the  linearized  solution,  but  as  nonlinearities 
become  important  the  perturbation  growth  accelerates  for  peaks  (a  >  0)  and  retards  for 
troughs  (a  <  0).  To  preserve  volume,  the  peaks  get  narrower  and  sharper  while  the 
troughs  broaden  and  flatten.  A  rough  estimate  for  the  duration  of  the  linear  behavior  is 
At  ~  \  In  1 1/ a0 1  - 

As  the  amplitude  gets  large,  the  growth  becomes  algebraic  in  the  time  remaining 
before  the  singularity  at  time  f*: 

ho  >  1  :  h~l/(t*—  t)  (4.3a) 

ho  <  1  :  h  « f*  —  t  (4.3b) 

which  shows  that  the  rapid  nonlinear  growth  occurs  over  a  time  scale  At  ~  1,  or,  dimen¬ 
sionally,  8/i/ Apgcli. 

The  catastrophic  growth  shown  by  the  inverse  relation  (4.3a)  between  peak  thickness 
and  remaining  time  is  strikingly  different  from  the  exponential  growth  of  small  ampli¬ 
tudes.  For  large  peaks  the  relevant  time  scale  is  inversely  proportional  to  the  current  peak 
thickness;  thus  large-amplitude  effects  drastically  enhance  the  growth  of  peaks. 

Another  effect  of  finite  amplitude  is  that  the  overall  strain  increases,  i.e.,  the  disturbed 
section  stretches  out  in  the  x  direction.  In  fact,  it  can  be  shown  (using  the  Schwarz 
inequality),  that  for  any  initial  perturbation  with  a  zero  mean,  the  length  of  the  perturbed 
section  increases  continually.  This  surprising  behavior  does  not  apply  in  general,  but  rather 
is  a  direct  consequence  of  the  assumption  that  the  force  on  the  ends  remains  constant. 

As  an  example,  if  the  initial  perturbation  is  sinusoidal  then: 


5o  =  1  +  bcosxo 


5(x0,t) 
x(xq, t ) 


2  +  6(1  +  e2t )  cos^o 
2  +  6(1  —  e2t )  cos^o 

— 6tanh(f)  sin(xo)  +  (sech2(f)  —  tanh(f))  xq 


2  tanh(f)sech(t) 

H -  =  arctan 

y  e~2t  —  6 2  cosh2(f) 


—  6cosh(f) 
+  6cosh(f) 


tan 


(4.4a) 

(4.4b) 


(4.4c) 


and  the  overall  strain  is  sech2(f)  —  tanh(t)  +  2  tanh(t)sech(t) / yj e~2t  —  b2  cosh2  ( t ) ,  which 
increases  monotonically  with  time.  The  layer  shapes  for  one  wavelength  of  this  solution 
at  various  times  are  shown  in  figure  3;  as  the  growth  becomes  nonlinear  the  wavelength 
stretches  out.  If  the  initial  amplitude  b  is  small,  the  solution  simplifies: 


1  +  |e2t  coscco 
1  —  |e2t  cosccq 


x 


A  -  (fa 


arctan 


tan(xo/2)  ]  —  x0 


(4.5a) 


(4.5b) 
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4.2.  CONSTANT  WAVELENGTH  DISTURBANCE 


For  the  second  case,  we  enforce  the  more  physically  reasonable  condition  that  the 
wavelength  of  the  disturbance  remain  constant  over  time.  (This  could  apply  either  to 
a  periodic  disturbance  of  infinite  extent,  or  to  a  layer  of  finite  extent  bounded  by  fixed 
shear-free  end  walls.)  Then  the  reduced  force  in  the  layer  must  vary  with  time  to  give 
zero  velocity  at  both  ends  of  one  wavelength.  Integrating  (3.11)  in  x  shows  that  fixed  ends 
require: 


(4.6) 


where  L  is  one  wavelength,  and  the  numerator  is  just  the  total  volume  of  fluid  in  the  layer 
over  that  wavelength,  and  so  remains  constant.  Thus  the  thinnest  parts  of  the  layer  have 
the  greatest  effect  on  the  current  equilibrium  thickness  de{t).  We  find  that  the  constant 
wavelength  requirement  causes  the  equilibrium  thickness  de(t)  to  decrease  in  over  time 
(the  proof  again  involves  the  Schwarz  inequality). 

Eliminating  de  from  (3.13)  (and  using  Sdx  =  Sodxo)  gives  a  single  partial  integro- 
differential  equation: 

DS  r2  Jo  dx o  .  . 

—  =  <r - p -  4.7) 

Dt  Jo  S0/S2dx o 

Again,  the  local  growth  depends  only  on  the  relation  between  the  local  thickness  and  an 
integral  property  of  the  entire  disturbance,  but  is  independent  of  the  immediate  neighbor¬ 
hood. 

For  large  amplitude,  the  main  effect  of  fixed  wavelength  is  to  slow  the  growth  of 
troughs;  in  fact,  the  thickness  6  is  prevented  from  reaching  zero.  In  contrast,  peak  growth 
is  only  slightly  accelerated.  When  peaks  become  large,  then  the  equilibrium  thickness 
de  becomes  negligible  in  comparison,  as  in  the  const  ant- de  case,  but  sooner  here  since 
de  decreases.  Thus  large  peaks  show  the  same  catastrophic  growth  to  infinite  thickness, 
given  by  (4.3a).  The  growth  of  large  peaks,  where  de  is  unimportant,  is  insensitive  to  end 
conditions,  and  our  previous  conclusion,  that  the  growth  of  peaks  is  dramatically  enhanced 
by  large- amplitude  effects,  applies  in  general. 

Qualitatively  this  case  is  very  similar  to  the  previous,  as  shown  by  the  profiles  in 
figure  4,  except  that  the  wavelength  remains  constant.  Figure  5  compares  the  growth  of 
the  disturbance  for  fixed  wavelength  with  that  for  constant  end  forces,  each  for  an  initial 
small  sinusoidal  disturbance. 


4.3.  LARGE-AMPLITUDE  BEHAVIOR 

Where  5  >  4,  we  can  approximate  (3.13)  by 


DS 

151 


=  S 2 


(4.8) 


(This  would  also  apply  if  the  layer  were  of  finite  extent,  surrounded  by  fluid  2,  so  de  —  0.) 
The  solution 

S(x0,t)  =  -  1  1  (4.9) 

<5o(®o) 
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shows  each  cross  section  reaching  infinite  thickness  in  a  finite  time  t  —  l/do‘,  this  is  the 
same  growth  as  (4.3a). 

For  this  case,  (3.11)  integrates  to: 

u(x,t)  =  -  6(£,t)d£  (4.10) 

Jo 

and  as  long  as  <5  remains  finite  everywhere,  we  can  express  this  in  Lagrangian  coordinates 
as: 

fXQ 

u(x0,t)  =  -  $o({,o)d£,0  (4.11) 

Jo 

Hence  each  cross  section  moves  at  a  uniform  speed  over  time.  These  results  apply  locally 
around  any  large  peak  where  cle  is  negligible. 

4.4.  AXISYMMETRIC  GROWTH 

How  does  the  three-dimensional  case  differ  from  the  two-dimensional  solutions  above? 
One  difference  is  that  the  reduced  force  is  a  tensor  that  varies  with  position,  rather  than  a 
uniform  scalar.  Consequently,  two  separate  portions  of  the  layer  having  the  same  thickness 
may  not  grow  at  the  same  rate.  Another  difference  is  the  possibility  of  (vertical)  vorticity 
and  shear  stresses  in  the  layer. 

Here,  we  examine  the  growth  of  an  axisymmetric  disturbance  without  swirl.  Then  the 
reduced  force  tensor  FtJ  is  diagonal  in  the  (r,  6)  coordinates,  without  shear,  and  has  radial 
and  circumferential  components 


Fn  =  62  +  5  (ur  +  \ u/r )  (4.12a) 

F22  —  +  <5  (lur  +  u/r)  (4.12b) 

where  u  is  the  radial  velocity.  The  radial  force  balance  (3.9)  can  be  expressed  as 

(Ui)r  +  P0)r  =  O  (4.13) 

The  axisymmetric  equations  (3.9,  3.10)  apparently  have  no  simple  closed-form  solu¬ 
tions,  except  when  the  layer  is  of  uniform  thickness.  In  this  special  case,  u  oc  r,  so  An  is 
uniform,  and  the  thickness  grows  according  to 

f  =  |(32:-Fn(0)  (4.14) 

Comparing  with  (3.13)  shows  that  the  axisymmetric  growth  of  a  flat  layer  is  qualitatively 
identical  to  (though  one-third  faster  than)  the  two-dimensional  case  for  comparable  force 
conditions  on  the  boundary.  This  special  case  may  be  indicative  of  the  dynamics  of  a 
smooth  peak  at  the  origin;  in  numerical  solutions  the  growth  of  large  peaks  approaches  the 
inverse  time  behavior  5  oc  (t*  —  t),  where  the  constant  of  proportionality  is  of  order  unity. 

We  have  solved  the  axisymmetric  equations  numerically  for  two  different  edge  condi¬ 
tions  and  a  variety  of  initial  conditions.  For  a  localized  disturbance,  at  the  edge  of  the 
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disturbed  disk  (ru)r  —  0  to  balance  the  force  in  the  surrounding  undisturbed  layer.  As  in 
the  two-dimensional  version,  the  disturbed  region  tends  to  spread  out  since  the  broadening 
of  the  troughs  overwhelms  the  narrowing  of  the  peaks.  For  a  disturbance  of  fixed  radial 
extent,  then  at  the  edge  u  =  0;  this  case  is  a  circular  approximation  of  one  hexagonal 
cell  of  a  periodic  array,  such  as  seen  in  the  experiments  of  Whitehead  and  Luther  (1975). 
Again,  the  spatial  constraint  limits  the  deepening  of  troughs,  while  slightly  accelerating 
the  growth  of  peaks.  In  figure  6,  we  show  profiles  for  this  case,  started  from  a  radial 
sinusoid.  This  shows  that  the  central  peak  grows  faster  than  the  ring-shaped  outer  peak, 
though  they  begin  with  identical  amplitudes. 

For  all  of  the  axisymmetric  cases  we  have  calculated,  the  qualitative  behavior  is  sim¬ 
ilar  to  the  two-dimensional  versions,  in  that  the  nonlinearity  leads  to  broad,  flat  troughs 
between  sharpening  peaks  whose  growth  accelerates  to  an  inverse-time  catastrophe.  One 
difference,  however,  is  that  the  growth  is  no  longer  merely  a  function  of  local  thickness, 
but  also  depends  on  position.  For  the  general  three-dimensional  case,  we  would  expect 
nonlinear  growth  to  yield  broad  depressions  separating  sharp  peaks  and  ridges  growing 
catastrophically  into  plumes  and  sheets,  although  there  could  be  additional  effects  due  to 
shear  and  swirl. 

5.  VERTICAL  SHEET  FORMATION 

Here  we  examine  what  happens  in  the  two-dimensional  case  when  a  thickness  max¬ 
imum  approaches  infinite  thickness  at  the  singular  time  f*.  We  will  show  how  the  peak 
becomes  a  vertical  sheet,  where  the  dense  viscous  fluid  drains  down.  The  small-slope  equa¬ 
tions  still  describe  the  behavior  in  most  of  the  layer,  except  in  an  asymptotically  narrow 
neighborhood  of  the  sheet,  because  the  sheet  has  only  an  asymptotically  small  effect  on 
the  stress  in  the  layer. 

As  a  peak  grows,  the  small-slope  assumption  must  break  down  locally  before  the  peak 
reaches  infinite  thickness.  Since  the  initial  disturbance  wavelength  is  long  compared  to 
the  layer  thickness,  then  when  the  peak  grows  larger  than  a  wavelength,  that  portion 
of  the  layer  around  the  peak  where  the  physical  slope  of  the  interface  is  0(1)  or  greater 
must  be  asymptotically  narrow  compared  to  the  whole  wavelength.  This  follows  from  mass 
conservation,  and  from  the  shape  near  the  peak  approaching  (as  we  will  show)  an  integrable 
negative  power  of  x  (where  x  =  0  at  the  peak).  Therefore  the  small-slope  approximation 
continues  to  apply  almost  everywhere  in  the  layer;  what  is  needed  is  a  description  of  the 
effect  of  the  peak  or  sheet  on  the  rest  of  the  layer. 

In  an  earlier  work  (Canright,  1987)  we  give  a  large-slope  analysis  appropriate  to 
the  region  around  a  peak  where  the  small-slope  approximation  no  longer  applies.  (The 
approach  used  there  is  essentially  that  of  Wilson,  1988.)  Assuming  the  interface  is  nearly 
vertical,  then  the  flow  is  extensional,  driven  by  negative  buoyancy  and  limited  by  normal 
viscous  stresses.  We  find  that  large-slope  effects  do  not  slow  down  the  growth,  they  only 
affect  the  details  of  the  peak  shape.  Indeed,  the  catastrophic  growth  described  by  (4.3a) 
still  applies  (except  for  a  numerical  coefficient  of  0(1)  depending  on  the  shape).  Physically, 
there  is  nothing  to  prevent  the  fluid  from  flowing  down,  and  so  the  peak  extends  to  become 
a  sheet.  Of  course,  at  some  point  the  extending  peak  either  will  reach  the  lower  boundary 
or  will  become  so  long  that  the  viscous  resistance  of  fluid  2  becomes  important.  In  the 
former  case,  a  pool  forms,  without  effect  on  the  upper  layer,  but  in  the  latter,  the  flow 
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driven  in  fluid  2  could  alter  the  dynamics  of  the  upper  layer. 

To  the  rest  of  the  layer  the  sheet  appears  as  an  isolated  singularity,  a  sink  of  fluid. 
The  horizontal  force  balance  must  still  apply  even  to  the  sheet,  and  so  the  reduced  force  in 
the  layer  F(t)  is  continuous  across  the  sheet.  In  general  F(t)  is  an  integral  property  of  the 
whole  layer,  for  example  (4.6)  for  a  fixed-wavelength  disturbance,  where  the  troughs  have 
a  greater  influence  than  the  peaks.  Because  the  large-slope  region  is  integrable  (as  it  must 
be,  since  no  new  mass  is  created)  and  asymptotically  thin,  its  effect  on  F(t)  is  negligible. 

As  an  example  of  the  formation  of  a  peak  singularity,  consider  what  happens  when  the 
initial  disturbance  is  a  small- amplitude  sinusoid.  (For  simplicity,  we  assume  the  constant- 
force  end  conditions,  but  since  F  has  little  effect  on  a  large  peak  the  results  apply  to  more 
general  end  conditions.)  The  previous  solution  (4.5)  shows  that  the  peak  becomes  singular 
as  be2t / 2  approaches  unity.  In  that  limit,  (4.5a, b)  become 


8  ~  cot2(a;o/2) 

(5.1a) 

x  ~  2  tan(xo/2)  —  Xq 

(5.1b) 

Then  near  the  peak  (x'o  <C  1) 

rp  3 

T  ~  xo 
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(5.1c) 

and  so 

5  m  (xq/2)~2  as  (3/2cc)— 2/3 

(5. Id) 

This  shows  that  at  the  singular  time  the  peak  becomes  proportional  to  an  integrable 
negative  power  of  x,  specifically  5  oc  x~ 2/3. 

In  fact,  the  same  power  of  x  results  from  any  smooth  initial  peak  that  locally  can 
be  fit  by  a  parabola.  As  the  singular  time  is  approached,  the  peak  can  be  described  by  a 
similarity  solution,  as  shown  in  Appendix  C.  The  general  similarity  solution  shows  that  a 
peak  of  the  more  general  form  50  ~  1  +  6(1  —  c|xo|n),  where  n  >  0,  gives  a  singularity  of 
the  form  S(x,  £*)  oc  where  m  =  n/(n  +  1).  Hence  any  initial  peak  becomes,  at  £*,  a 

singularity  proportional  to  an  integrable  negative  power  of  x. 

To  follow  how  fluid  flows  into  the  sheet,  the  large-amplitude  equations  (4.9,  4.10)  apply 
near  the  sheet  (but  outside  the  large-slope  region).  After  t*,  there  is  a  singularity  (the 
sheet)  at  the  origin,  and  the  layer  must  move  in  such  a  way  that  each  fluid  cross-section 
reaches  the  singularity  at  the  same  time  that  it  reaches  infinite  thickness.  This  determines 
the  strength  of  the  mass  sink  over  time. 

Consider  the  fate  of  the  inverse  power  of  x  singularity  that  forms  at  f*.  This  is  most 
simply  described  if  we  take  the  state  of  the  layer  at  f*  to  be  the  new  reference  state, 
and  relabel  each  fluid  cross  section  x0  by  its  new  reference  position  cc*  =  x(xq ,  t*),  with 
reference  thickness  h*(a;*)  =  8(x,t *),  so 

at  r  =  0  :  5*(a;*)  =  S(x  —  x*, t  —  f*)  =  bx~m  ,  0<m<l  (5.2a) 

where  r  =  t  —  f*  and  b  is  some  positive  amplitude.  Then  the  fluid  cross  section  which 
arrives  at  the  singularity  at  time  r  is  that  which  becomes  infinitely  thick  at  that  time;  we 
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call  this  section  xs(r)  :  x(x *  =  xs(r),r)  =  0.  Using  (4.9)  gives 


xa(r)  =  ( br 

(5.2b) 

S(x*,T)  =  b/(x?-x?) 

(5.2c) 

and  so 

x(x*,  T )  =  [  \d£*  =  (xm  xa)  s  (xl  m  x]  m) 

Jxa  o(s*>  t)  1-m 

(5. 2d) 

Away  from  the  sheet  (x*  3>  xs),  the  profile  is  still  the  starting  profile  from  r  = 

* 

ll 

o 

8  ~  bx  m) 

.  However,  very  close  to  the  sheet: 

for  e  =  {x*/xs  —  1)  <C  1 

(5.2e) 

x  ~  mxse2 / 2 

(5.2f) 

8/b  a;  1  /(mex™)  ~  j \/ 2mx 

(5.2g) 

This  shows  that  very  near  the  sheet,  the  singularity  goes  like  l/y^,  with  a  scale 
that  varies  in  time.  (This  local  x  dependence  is  actually  independent  of  the  starting 
conditions,  as  shown  in  Canright,  1987,  App.  C.)  From  (5.2g)  it  is  clear  that  whether 
the  thickness  around  the  sheet  grows  or  shrinks  is  determined  by  whether  m  <  1/2  or 
m  >  1/2,  respectively. 

The  special  case  m  —  1/2  gives  a  steady  solution.  For  m  >  1/2,  the  fluid  drains 
away  down  the  sheet  faster  than  it  comes  in  from  the  sides,  and  the  square-root  singularity 
diminishes  with  time  as  it  spreads  out,  to  match  onto  the  nearly  undisturbed  profile  x~m. 
This  would  be  the  eventual  fate  of  an  initially  ( t  =  0)  smooth  maximum,  which  gives 
m  —  2/3.  Conversely,  if  m  <  1/2,  the  square-root  singularity  grows  as  it  spreads,  fed  from 
the  sides  faster  than  it  can  drain  fluid  away.  (To  get  m  <  1/2  would  require  a  cusp-like 
initial  [t  —  0]  maximum,  which  may  not  be  physically  realistic.)  This  solution  (5.2)  is 
again  a  particular  case  of  the  general  large-amplitude  similarity  solution  of  appendix  C. 

(For  the  axisymmetric  case,  there  is  no  similarity  solution  that  describes  how  a  finite 
peak  grows  into  a  plume,  but  we  speculate  that  the  same  qualitative  behavior  applies.  A 
steady  axisymmetric  plume  has  the  shape  8  oc  1/r.) 

With  the  above  description  of  how  a  sheet  first  forms  and  how  it  behaves  afterward, 
the  small-slope  equations  can  be  used  to  follow  the  development  of  the  instability  from 
initial  conditions  through  rapid  large-amplitude  growth  all  the  way  to  the  draining  away 
of  the  fluid  down  the  sheets.  The  results  will  be  inaccurate  wherever  the  physical  slope  of 
the  interface  is  not  small,  but  such  regions  comprise  only  a  small  fraction  of  the  domain 
and  have  little  effect  on  the  dynamics  of  the  rest  of  the  layer.  The  only  assumption 
is  that  a  sheet  does  not  exert  any  net  horizontal  force  on  the  surrounding  layer.  This 
assumption  may  break  down  if  the  length  of  a  plume  becomes  so  much  greater  than  the 
initial  wavelength  that  the  flow  it  drives  in  the  lower  fluid  becomes  dynamically  significant. 

6.  CONCLUSIONS 

The  central  concern  of  this  work  is  the  nonlinear  interactions  between  buoyant  forces 
and  normal  viscous  stresses  that  occur  in  a  buoyantly  unstable  viscous  layer  under  a  shear- 
free  horizontal  boundary  and  over  a  much  less  viscous  fluid.  After  the  initially  uniform 
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thickness  of  the  layer  has  been  perturbed  slightly,  the  early  growth  of  the  perturbation 
is  exponential;  the  perturbation  keeps  its  shape  while  it  grows.  But  when  the  nonlinear 
effects  of  finite  amplitude  become  important,  the  thicker  parts  of  the  layer  thicken  more 
rapidly  while  the  thinner  parts  thin  more  slowly,  giving  in  general  sharp  peaks  with  broad, 
flat  troughs  in  between,  over  a  timescale  of  8/x/ Ap  gd\.  The  accelerating  growth  of  peaks 
leads  to  infinite  thickness  at  some  time  t*,  and  the  final  catastrophic  growth  of  the  peak 
thickness  5  is  algebraic:  5  ~  fj,/ Apg(t*—t)  for  the  two-dimensional  case.  The  axisymmetric 
case  shows  essentially  the  same  behavior.  Similar  catastrophic  growth  is  also  predicted  for 
a  power-law  fluid,  though  the  power  of  (t*  —t)  and  the  coefficient  are  different.  This  shows 
how  large- amplitude  growth  is  fundamentally  different  from  small-amplitude  growth;  large- 
amplitude  effects  dramatically  enhance  the  growth  of  peaks. 

The  small-slope  equations  continue  to  apply  to  the  layer  even  after  the  formation  of 
downwelling  sheets,  except  in  an  asymptotically  narrow  neighborhood  around  each  sheet. 
This  is  possible  because  the  sheets  do  not  change  the  horizontal  force  balance  (unless  the 
flow  they  drive  in  the  lower  fluid  becomes  dynamically  significant).  Applying  the  equations 
up  to  the  singular  time  shows  that  at  first  the  sheet  should  have  the  local  shape  5  oc  |x|-2/3, 
but  that  afterwards,  as  the  sheet  drains  the  layer,  the  sheet  changes  shape  to  6  oc  \x\~1^2. 
This  behavior  is  clarified  by  a  family  of  similarity  solutions,  appropriate  where  5  is  large. 

This  analysis  depends  on  two  key  assumptions:  small  interfacial  slope  and  negligible 
shear  stress  from  fluid  2.  When  peaks  become  large  enough,  the  slope  becomes  large, 
and  the  approximate  equations  become  invalid.  However,  as  long  as  the  external  shear 
is  negligible,  the  growth  of  the  disturbance  at  large  slopes  is  essentially  the  same  as  that 
predicted  here  (Canright,  1987,  or  Wilson,  1988),  with  catastrophic  inverse-time  peak 
growth;  while  the  details  of  the  peak  shape  are  different,  this  only  modifies  the  prediction 
of  f*  by  an  0(1)  numerical  factor.  Even  so,  when  the  descending  sheet  becomes  long 
enough,  the  flow  driven  in  fluid  2  will  result  in  significant  shear  stress  on  the  interface, 
retarding  the  growth  and  invalidating  the  approximate  equations.  (It  is  also  possible, 
depending  on  the  parameters,  that  inertia  could  become  important  for  rapidly  growing 
peaks,  or  that  surface  tension  could  become  significant  at  the  highly  curved  peaks.)  So 
this  analysis  is  not  appropriate  to  describe  what  happens  at  the  tip  of  a  long  descending 
sheet;  e.g.,  the  tip  might  widen  to  a  bulb,  as  seen  in  the  experimental  plumes  of  Whitehead 
and  Luther  (1975)  and  the  numerical  sheets  of  Newhouse  and  Pozrikidis  (1990)  (though 
both  of  those  works  have  a  no-slip  surface,  a  significant  difference).  Rather,  the  present 
analysis  describes  the  transition  from  exponential  to  catastrophic  growth  of  peaks  due  to 
finite  amplitude,  and  the  subsequent  adjustment  of  the  rest  of  the  layer. 

The  question  of  whether  plumes  or  sheets  are  more  likely  to  develop  is  beyond  the 
scope  of  this  work.  Indeed,  the  dynamic  balance  considered  is  independent  of  wavelength, 
and  so  presumably  independent  of  planform  as  well,  at  least  for  Newtonian  fluids.  However, 
we  speculate  that  the  lack  of  downwelling  plumes  in  the  mantle  is  due  to  rheological  effects. 
Assuming  the  lithosphere  weakens  with  deformation,  the  two-dimensional  motion  leading 
to  sheets  allows  the  deformation  to  be  concentrated  at  the  growing  peak  while  the  layer 
on  both  sides  moves  rigidly  inward;  the  axisymmetric  flow  needed  for  plumes  necessitates 
significant  deformation  throughout  the  entire  layer. 

To  apply  these  results  to  the  oceanic  lithosphere,  an  order-of-magnitude  estimate 
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of  the  time  scale  for  the  formation  of  a  new  subduction  zone  (starting  from  a  finite- 
amplitude  disturbance,  say  due  to  mantle  convection)  by  this  process  is  about  0.3  Gyr, 
or  one-tenth  the  age  of  the  earth.  (This  assumes  /ii  ~  1025  poises  [Walcott,  1973],  d\  ~ 
100  km,  and  A p  ~  O.lg/cm3  from  a  temperature  difference  AT  ~  1000  K  with  a  thermal 
expansion  coefficient  a  ~  3  x  10-5  /K  from  a  base  density  of  p  ~  3.3g/cm3.)  This  figure 
is  for  illustrative  purposes  only;  the  uncertainty  in  the  appropriate  viscosity  is  orders  of 
magnitude.  However,  if  the  surface  viscosity  of  the  early  earth  were  orders  of  magnitude 
smaller  than  today,  as  has  been  suggested  in  hot-earth  models  (Davies,  1990),  then  the 
time  scale  would  be  correspondingly  reduced;  one  possible  interpretation  would  be  that 
modem  subduction  zones  may  have  their  origins  early  in  the  earth’s  history  when  the 
surface  layer  was  sufficiently  deformable. 


APPENDICES 

Appendix  A.  LINEARIZED  SOLUTION 

Here  we  calculate  the  small-amplitude  growth  rates  for  the  non-iuertial  problem  with 
arbitrary  wavelength,  depths,  and  viscosities.  (This  is  an  extension  of  the  analysis  given 
by  Whitehead  and  Luther,  1975,  to  include  the  effects  of  finite  depth  in  fluid  2.)  A  more 
convenient  reduced  pressure  p  in  each  fluid  is  defined  in  terms  of  the  total  pressure  P  by 


Pi=Pi+Pigz  (A.  la) 

P-2  =  P2  +  p-2gz  +  Apgch  (A.  lb) 


(or  from  the  other  reduced  pressure  p  by  pi  =  pi  - \-Apgz ,  p2  =  p2+^pgdi)  so  in  equilibrium 
p  is  zero.  Assuming  the  interface  disturbance  varies  only  in  the  x  direction  and  has  small 
amplitude  and  slope,  the  resulting  2-D  problem  can  be  linearized  by  applying  the  interface 
conditions  at  z  —  d\ .  Then  the  interface  conditions  (2.4b-d)  become: 


at  z  —  d\  :  [it]  =  \p{uz  +  wx)]  =  0  (A. 2a) 

[— P  +  2nwz]  =  Apg(6  -  d])  (A. 2b) 

5t  —  w  (A. 2c) 


2 

where  again  []|“  indicates  the  jump  in  value  from  fluid  1  to  fluid  2.  This  problem  is 
separable  and  has  a  simple  analytical  solution,  assuming  the  perturbation  is  sinusoidal. 

The  solution  is  given  below,  where  k  is  the  wave  number,  a(t)  is  the  dimensionless 
amplitude,  Z  =  z  —  (di  +  c^)  is  the  coordinate  in  fluid  2,  a  =  P2/ Pi  is  the  viscosity  ratio, 
and  in  the  coefficients  A,  B ,  E.  F.  and  their  common  denominator  D  these  abbreviations 
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are  used  :  k  =  2kdi,  K  =  2/cd2,  c  =  cosh(/c),  (7  =  cosh(JF),  s  =  sinli(fc),  S'  =  sinh(Ji). 


5  —  d\  [1  +  a(t )  cos(for)]  (A.3a-j) 

d/i  =  (Apgdi/nik2)a(t)/2sm(kx)[Asmh(kz)  +  Bkzcosh(kz)] 

\I/ 2  =  {Apgdi/ pik2)a(t) /2  sm(kx)[E  smh(kZ)  +  FkZ  cosh(kZ)] 
pi  =  (Apgdi){—a(t)B  cos(kx)  cosh(kz)} 
p2  =  (  A/9(7cZi){—  a(t)aF  cos(kx)  cosh(fcZ’)} 

A  =  {-1/D){[2(S  -  K)  +  ak{C  -  1)]  sinh(fcdi) 

+  [k(S  —  K )  +  a(kK  +  2 (C  —  1))]  cosh(kdi)} 

B  =  (2 /D){(S  -  I<  +  aK)  sinh(fcdi)  +  a(C  -  1)  cosh(fcdi)} 

E  =  (1/D){[2a(s  -  k)  +  K(c  -  1)]  sinh(fcd2) 

+  [aK (s  —  k)  +  Kk  +  2(c  —  1)]  cosh(/cd2)} 

F  =  (— 2/D){[a(s  —  k)  +  k]  shih(kd2)  +  (c  —  1)  cosh(/cd2)} 

D  =  (S  —  K)(s  +  k)  +  2 a(Cc  -  1  +  Kk)  +  a2(S  +  K)(s  -  ic) 


Then  from  St  =  w(z  —  dp)  we  get  the  growth  rate: 

a(t)  —  a(0)ecrt 
a  =  [Apgdx/pPjd 

1  (S -K)(c-l)  +  a(s-k)(C -1) 

k{S-  K)(s  +  k)  +  2 a(Cc  -  1  +  Kk)  +  a2(S  +  K)(s  -  k) 


(A. 4a) 
(A. 4b) 

(A. 4c) 


The  symmetry  of  the  problem  is  apparent  in  the  solution.  Thus,  without  loss  of  generality, 
assume  that  fluid  2  is  the  deeper  layer:  K  >  k. 

This  solution  is  governed  by  three  dimensionless  parameters:  the  non-dimensional 
wavenumber  k  (or  e),  the  depth  ratio  j3  =  d2/di  =  k/k,  and  the  viscosity  ratio  a  =  p2/pi. 
Note  that  d  is  a  monotonically  decreasing  function  of  ct;  if  we  increase  /r2  while  pi  stays 
constant  the  growth  rate  can  only  decrease.  In  the  limit  f3  — >  oo  then  (A. 4c)  reduces  to: 


,  =  i  (c_i)  +  r(s-fc)  _  (A8) 

k  (s  +  k )  +  2rc  +  r2(s  —  k ) 

which  is  just  the  result  of  Whitehead  and  Luther  (1975).  Figure  7a  shows  the  effects  of 
finite  depth  for  aCl;  figure  7b  shows  the  effects  of  viscosity  ratio  on  d(k)  for  j3  —  oo. 

When  /f  >  1  there  are  well  defined  regimes  of  growth  where  different  force  balances 
are  dominant.  These  are  shown  schematically,  with  the  corresponding  growth  rates,  in 
figure  8.  In  discussing  the  various  growth  mechanisms  below,  it  should  be  kept  in  mind 
that  the  same  mechanisms  continue  to  apply  as  long  as  the  slope  of  the  interface  remains 
small,  which  for  long  waves  includes  large-amplitude  growth. 

For  sufficiently  short  waves  (k  3>  min(l, max(/3-1, ck-1/3)))  the  disturbance  sees  nei¬ 
ther  boundary  and  a  [k(  1  -Feu)]-1,  so  if  one  viscosity  is  much  larger,  that  one  limits  the 
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growth.  The  dominant  mechanism  here  is  that  the  vertical  motion  of  the  interface  is  re¬ 
sisted  by  normal  viscous  stresses  in  the  more  viscous  fluid,  and  the  growth  rate  diminishes 
with  decreasing  wavelength. 

At  the  other  extreme,  for  sufficiently  long  waves  (k  -C  min(/3-1,  y/ft/a,  \ /a/?))  the 
boundaries  confine  the  flow  to  be  mainly  horizontal,  limited  by  shear  stresses  at  the  inter¬ 
face.  Then  a  — >  k2(  1  +  and  the  controlling  viscosity  depends  on  whether  a  >  (5 

or  not. 

Between  these  extremes,  the  waves  are  long  compared  to  layer  1  so  the  motion  of 
the  interface  is  primarily  horizontal,  and  there  are  four  different  regimes.  In  one  (a-1  -C 
k  <  a-1/3,  l«a<  /33),  fluid  2  is  relatively  immobile  and  the  growth  is  limited  by  the 
shear  across  layer  1,  giving  the  same  growth  rate  as  the  previous  case,  i.e. ,  k2  j  12.  (This 
possibility  was  apparently  overlooked  by  Whitehead  and  Luther.)  In  another  (/?_1  -C  k  <C 
min(a,a-1)),  the  slight  resistance  of  fluid  2  (the  rate-controlling  viscosity  is  /j, 2)  gives  a 
small  shear  gradient  across  layer  1,  which  over  the  long  wavelength  is  sufficient  to  balance 
the  buoyancy,  giving  a  — >  k/ 4a. 

In  the  other  two  regimes,  the  less  viscous  fluid  is  effectively  passive,  and  the  primarily 
horizontal  motion  of  the  more  viscous  layer  is  limited  by  normal  viscous  stresses.  The 
resulting  growth  rate  is  nearly  independent  of  wavelength,  and  includes  the  maximum 
growth  rate  possible  for  a  given  viscosity  contrast  a  where  this  behavior  occurs.  (For 
other  a,  i.e.,  1  <C  a  <  /l3,  we  expect  the  fastest  growth  at  the  crossover  between  short- 
and  long-wave  behavior.)  When  fluid  2  is  much  more  viscous  (a  3>  /33),  this  regime 
(yJ/3/a  <C  k  -C  /3_1)  gives  a  — >  /3/4a,  while  for  fluid  1  more  viscous  (a  -C  1)  this  regime 
(max(a,  y/af3)  <C  1)  gives  <7  — >  1/4. 

In  the  last  case,  examining  the  broad  peak  more  closely  shows  that,  for  /3~5  <a<l, 
a  ~  ( 1/4) ( 1  —  (a/k  +  k4 /720)).  This  broad  maximum  peaks  at  kmax  ~  (180a)1/5  =  2.8a1/5 
and  amax  ~  ( 1/ 4) ( 1  —  0.44a4/5).  As  an  indication  of  the  flatness  of  this  peak,  for  the  range 
a4/5  <  k  <  5.2a1/20,  a  >  (1/4) (1  —  a1/5).  When  a  • <  /3~5 ,  the  finite  depth  modifies 
the  maximum  growth  rate  giving  a  ~  (1/4)(1  —  (3 a//3k2  +  /c4/ 720))  with  a  broad  peak  at 
kmax  ~  3.2(a/ /3)1/6  and  amax  ~  (1/4/1  -  0.44(a//)2/3). 

The  present  work  is  only  concerned  with  the  case  of  a  thin,  viscous  layer  over  a  less 
viscous,  deep  layer,  so  a  C  1  and  j3  3>  1.  We  further  restrict  our  consideration  to  the 
mechanism  giving  the  fastest  growth,  i.e.,  the  last  regime  considered  above,  where  a  ~  1/4. 
The  lowest-order  finite- amplitude  analysis  (section  3)  therefore  predicts  that  the  growth  is 
independent  of  wavelength. 

Appendix  B.  POWER-LAW  FLUID 

The  long-wave  analysis  of  section  3  is  not  limited  to  Newtonian  fluids.  The  important 
point  is  that  both  surfaces  of  the  layer  are  shear-free,  so  throughout  the  layer  shear  stresses 
are  O(e)  smaller  than  normal  stresses,  and  the  latter  are  independent  of  2  to  0(e2).  (This 
assumes  the  lower  fluid  is  passive  and  effectively  hydrostatic;  the  restriction  this  applies  to 
wavelength  will  depend  on  the  rheology  and  depth  of  fluid  2.)  Hence  the  horizontal  force 
balance  (3.5,  3.6)  applies  for  any  rheology. 

Here  we  consider  a  layer  of  fluid  with  a  power-law  rheology,  subject  to  a  disturbance 
independent  of  y.  Then  the  constitutive  relation  can  be  written  in  the  form: 

Txx  =  TZz  =  %  Hr  Sgn(lta;) \ux  \  /  (-^-1) 
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where  fir  has  appropriate  units  and  m  >  1  for  a  fluid  that  weakens  with  increasing  strain 
rate.  Eliminating  rxx  and  u  as  before  yields  a  single  Lagrangian  equation  describing  the 
evolution  of  the  thickness  of  the  layer  as  we  follow  a  material  cross  section,  in  dimensionless 
form: 


DS 

Dt 


sgn(5  -  de(t))S 


5  -d2e{t)  | 

5 


(B.2) 


using  the  same  scales  as  before, 
special  case  of  a  Newtonian  fluid 
we  will  assume  that  m  is  an  odd 
For  an  infinitely  long  layer 

de  =  1) 


except  for  the  new  time  scale  (8/ir/ Apgd\)m .  For  the 
(m  =  1)  this  reduces  to  (3.13).  To  simplify  the  notation, 
integer. 

with  a  localized  disturbance  (i.e.,  constant  end  forces: 


D5  _  ( 5 2  -  l)m 
Dt  ~  hm_1 


(B.3) 


This  can  be  integrated  by  parts.  For  example,  if  m  =  3: 


r(5) 


5  5  1 

4(52  _  !)2  +  8(£2  _  1)  +  16  111 


5-1 

J+l 


(B.4) 


where  r(5)  =  t*  —  t  is  the  time  remaining  before  the  thickness  of  the  fluid  cross  section 
goes  to  oo  or  0. 

While  the  amplitude  a  =  d  —  1  remains  small,  the  approximate  solution  is: 


ao 

[1  -  (m  —  1/(m-I) 


(B.5) 


where  ao  is  the  initial  amplitude.  (This  algebraic  growth  is  quite  different  from  the  small- 
amplitude  exponential  growth  for  Newtonian  rheology.)  This  slow  growth  lasts  for  a  period 
of  roughly  At  ~  l/[(m  —  l)2maon_1].  For  large  amplitudes,  the  growth  again  becomes 
algebraic  in  the  time  remaining  before  the  singularity  is  reached  at  t  =  £*: 

50  >  1  :  5  «  l/[m(f*  —  t)]1/™  (B.6a) 

50  <  1  :  5  «  [m(U  -t)]1/m  (B.6b) 


This  catastrophic  growth  occurs  over  a  time  scale  At  ~  1/m  (or  in  dimensional  terms 
(8/xr/  Apgh)m /m). 

As  in  the  Newtonian  case,  a  disturbed  region  under  constant  end  forces  will  tend  to 
stretch  out  in  the  x  direction.  To  keep  the  wavelength  constant,  the  end  forces  must  vary 
in  time  to  give: 


[L 

\5-dim 

Jo 

5 

dx 


0 


(B.7) 


Layer  profiles  calculated  using  this  condition  are  shown  in  figure  9  for  m  =  3  and 
m  —  7,  from  an  initial  sinusoidal  disturbance.  Profiles  for  the  constant-end-force  condition 
(not  shown)  are  qualitatively  similar.  The  effect  of  increasing  m  is  seen  to  be  concentration 
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of  the  deformation  into  narrow  regions  at  the  centers  of  peaks  and  troughs,  while  elsewhere 
the  profile  becomes  linear  in  x. 

Figure  10  shows  how  the  initial  slow  growth  suddenly  becomes  catastrophic  after 
a  certain  threshold  has  been  reached  (more  so  for  higher  m).  For  the  constant-force 
condition,  a  trough  reaches  this  threshold  and  necks  off  much  sooner  than  a  peak  of  equal 
initial  amplitude  blows  up.  The  constant-wavelength  condition,  however,  causes  peaks  to 
grow  sooner  than  troughs. 

Regardless  of  the  end  conditions,  when  peaks  get  large  (compared  to  the  current  equi¬ 
librium  thickness),  the  large-amplitude  effects  still  produce  catastrophic  growth  algebraic 
in  the  time  remaining  (S  oc  l/(£*  —  t)1//m),  giving  infinite  thickness  in  finite  time,  as  for  a 
Newtonian  fluid. 

Appendix  C.  LARGE- AMPLITUDE  SIMILARITY  SOLUTION 

The  equations  appropriate  to  large-amplitude  disturbances  (for  the  two-dimensional 
case)  admit  a  rich  family  of  similarity  solutions,  which  illustrate  a  variety  of  behaviors. 
While  such  solutions  demand  particular  initial  conditions,  nonetheless  these  solutions  can 
be  interpreted  as  good  local  approximations  for  situations  arising  from  arbitrary  initial 
conditions.  Two  cases  are  of  particular  relevance  to  sheet  formation:  one  describes  how 
a  smooth  finite  peak  evolves  to  an  infinite  singularity,  the  other  describes  how  that  first 
singularity  changes  shape  as  the  sheet  evolves. 

For  large  amplitudes  d(  is  negligible,  and  the  Eulerian  equations  are: 

5t  +  ( u5)x  —  0 

8  +  ux  =  0 

We  assume  a  similarity  solution  of  the  following  form: 

x  =  a(r)f 

5  =  A(r)m  (C.2) 

u  =  a{r)A{r)g(Q 

where  r  =  t  — f*  is  the  time  relative  to  the  singular  time  f*.  Substition  shows  that  similarity 
requires  both  A') A2  and  a'/(ciA )  to  be  constants. 

Choosing  the  scale  of  A  (provided  A'  ^  0)  gives 

A(t)  =  1/r  (C.3a) 

Of  course,  <5  must  be  non-negative,  so  solutions  with  f  >  0  will  apply  only  for  r  >  0, 
when  A  is  decreasing,  and  those  with  /  <  0  will  apply  for  r  <  0,  when  —A  is  increasing 
catastrophically  to  r  =  0.  The  second  constant,  say  A,  implies: 

a(r)  =  |r|A  (C.3b) 

For  r  >  0,  positive  A  gives  a  profile  that  spreads  out  in  the  x  direction,  while  for  negative 
A  the  profile  contracts;  the  converse  is  true  for  r  <  0. 


(C.la) 

(C.lb) 
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The  system  (C.l)  becomes 

“/ -A£/' +  (/<?)' =  0  (C.4a) 

f  +  9'  =  0  (C.4b) 

The  solution  involves  two  arbitrary  constants  C  and  D  and  a  change  of  variables  (as  long 
as  A  /  0.  ij,  corresponding  to  a  different  inertial  reference  frame: 

C  =  £  -  C/{ A2  -  A) 

g(C)  =  g(t)  -  C/{\  -  1)  (C.5) 

Then  the  solution  is  given  implicitly  by 

C  —  q  +  sgn(C  —  q)D\q\k  (C.6a) 

/  =  — !/[l  +  sgn(g(C  -  q))kD\q\k~l]  (C.6b) 


where  k  =  A/(A  —  1)  and  D  >  0.  This  is  the  general  solution;  a  phase-plane  analysis  verifies 
that  this  gives  all  solutions  (for  A(t)  not  constant,  A  ^  0, 1)  except  the  trivial  solutions 
/  =  0  and  /  =  —1.  This  solution  describes  a  variety  of  behaviors,  depending  primarily 
on  A  and  on  which  branch  of  the  solution  is  chosen.  The  constants  C  and  D  affect  the  x 
origin  and  scale,  respectively. 

For  example,  consider  A  >  1  (so  k  >  1),  C  =  0,  D  —  1,  and  the  branch  where 
sgn(£  -  g)  =  sgn(£)  =  sgn (g): 

£  =  g  +  sgn(g)\g\k  (C.7a) 

f  =  -!/[! +  k\g\k-1}  (C.7b) 

Since  /  <  0,  this  solution  only  applies  for  r  <  0;  it  describes  the  growth  of  a  finite  peak  up 
to  the  singular  time  when  6  becomes  infinite.  The  asymptotics  in  £  reveal  the  behavior: 

[l-fc|r|-fc|x|fc-1]/|r[ 

—x/\t\  (C.7c) 

|  x\-xlxlk 

— sgn(x)|x|1/fc  (C.7d) 

The  first  (C.7c)  applies  for  x  near  the  origin  (|x|  <C  |t|a).  The  peak  is  of  fairly  general 
shape,  but  to  be  analytic  in  x,  k  must  be  an  odd  integer.  The  second  (C.7d)  applies  far 
from  the  origin  early  on  (r  0),  but  applies  ever  nearer  until  at  the  singular  time,  it 
applies  for  all  x.  The  asymptotic  shape  (C.7d)  is  independent  of  time;  as  the  peak  grows  it 
fills  in  the  profile  of  the  integrable  negative  power  of  x.  For  a  smooth  peak,  k  —  3,  A  =  8/2 
and  at  the  singular  time  5  oc  |a;|-2/3  (for  any  D). 

As  a  second  example,  A  >  1,  C  =  sgn(£),  D  —  (A  —  l)fc/A,  and  the  branch  where 
sgn(0  =  -sgn(fif): 

i  =  g-  sgn(g)A-1[((A  -  l)\g\  +  l)fc  -  1] 

/  —  1/ [(( A  —  l)\g\  +  l)fc_1  —  1] 


as  |^|  — >  0  :  5 

u 

as  |£|  — >  oo  :  5 

u 
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(C.8a) 

(C.8b) 


which  describes  a  sheet,  symmetric  in  x,  for  r  >  0.  Asymptotically: 

5— >1/[t1-a/2v/2^[] 

u  — >  — sgn(a;)-\/2|:r|/T1_A/2  (C.8c) 

&  -»•  (AM)“1/A 

u  — >  —  sgn(x)(A|a:|)1/fc/(A  —  1)  (C.8d) 

At  the  singular  time  (r  =  0),  the  time-independent  profile  (C.8d)  applies  everywhere;  this 
example  shows  what  happens  to  an  initial  profile  proportional  to  an  integrable  negative 
power  of  x.  (By  a  different  choice  of  C  and  D,  the  profile  (C.8d)  could  be  made  to  match 
(C.7d)  exactly.)  For  r  >  0,  (C.8c)  applies  near  the  origin;  regardless  of  the  starting  x 
dependence,  the  shape  around  the  singularity  (sheet)  becomes  proportional  to  l/yf\x[. 
If  A  <  2  then  the  square-root  singularity  and  the  corresponding  strength  of  the  sink  at 
the  origin  decay  with  time,  as  the  layer  drains  away.  (This  would  be  true  for  an  initially 
smooth  profile  like  (C.7c)  with  A  =  3/2.)  Conversely,  for  A  >  2,  the  singularity  grows,  as 
fluid  comes  in  from  the  sides  faster  than  it  can  be  disposed  of.  The  special  case  A  =  2 
gives  a  steady  sheet.  There  are  also  related  solutions  for  A  =  0,  0  <  A  <  1,  and  A  =  1 , 
which  have  the  same  asymptotic  shapes  (C.8c,d),  except  that,  for  |f  |  — ►  oo,  the  forms  for 
u  are  different,  and  for  A  =  0,  5  — >  e-^ / 1  as  | ^ |  — ^  oo . 

Briefly,  the  other  behaviors  governed  by  the  similarity  solution  are  as  follows.  For 
t  <  0,  i.e.,  profiles  growing  to  the  singular  time,  there  are  four:  (i)  a  zero  minimum  of 
shape  \x\k ,  k  >  0,  locally  time-independent,  that  far  away  levels  off  to  approach  l/|r|;  (ii) 
a  profile  that  approaches  zero  as  x  — >  —  oo  and  l/|r|  as  x  — >  +oo;  (iii)  a  symmetric  finite 
minimum  flanked  by  sheets  (which  could  be  extended  periodically);  and  (iv)  a  sheet  whose 
sides  level  off  to  approach  l/|r|,  rather  than  zero.  For  r  >  0,  where  the  profile  starts  at 
the  singular  time  and  diminishes  thereafter,  there  is  only  one  other  case:  a  zero  minimum 
flanked  by  sheets  (which  could  extend  periodically).  The  steady-shape  case  of  A(t)  =  1,  so 
a(r)  —  eAr,  gives  either  an  isolated  sheet,  or  a  zero  minimum  where  |5X|  — >  oo  surrounded 
by  sheets  (possibly  periodic).  Note  that  solutions  with  the  same  A{r)  and  A  but  different 
C  and  D  can  be  “cut  and  pasted”  together  to  give  other  similarity  profiles;  if  g  and  g' 
are  continuous  this  will  produce  reasonable  results.  All  cases  with  sheets  show  the  same 
1/ \f\x\  shape. 
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FIGURE  CAPTIONS 


Figure  1.  Problem  description.  Between  two  horizontal,  shear-free  planes,  a  layer  of  fluid  1 
of  uniform  thickness  d\  is  initially  in  unstable  equilibrium  over  a  layer  of  fluid  2  of  thickness 
d2.  The  upper  layer  is  denser  ip\  >  p2),  much  thinner  (d\  <C  d2),  and  much  more  viscous 
(pi  3>  y2).  After  the  interface  is  perturbed,  its  position  is  given  by  5(x,y,t). 

Figure  2.  Balance  of  horizontal  forces,  reduced  by  the  hydrostatic  pressure  gradient  of 
fluid  2,  on  a  differential  column  of  fluid  1.  Because  there  is  no  shear  above  or  below  and 
no  (reduced)  pressure  below,  the  reduced  horizontal  force  tensor  F^j  acting  on  the  vertical 
surfaces  of  the  column  has  zero  divergence. 

Figure  3.  Interface  profiles  5(x)  at  various  times  ( t  —  2.5,  3.0,  3.5,  3.7)  following  the 
evolution  of  one  wavelength  of  an  initial  sinusoidal  perturbation  of  amplitude  10-3,  for 
the  constant-end-force  condition  (F  =  1).  (Here  the  peak  reaches  infinite  thickness  at 
f*  =  3.80.)  As  nonlinear  effects  become  important,  the  peak  sharpens  while  the  troughs 
flatten  and  widen;  also  the  wavelength  increases  over  time. 

Figure  4.  Interface  profiles  for  the  constant-wavelength  end  condition  (L  =  1).  Similar 
to  Figure  4,  except  that  the  troughs  become  thin  more  slowly,  the  peaks  reach  infinite 
thickness  slightly  sooner  (at  t*  =  3.74),  and  the  wavelength  remains  constant. 

Figure  5.  Growth  of  the  disturbance  amplitude  over  time  for  an  initial  perturbation 
amplitude  of  10-3.  Both  the  peak  amplitude  amaa;  =  8max  —  1  and  the  trough  am¬ 
plitude  amin  =  1  —  hmin  are  shown:  constant-wavelength  condition  (solid),  constant 
end- force  (dashed),  and  linearized  (exponential)  growth  (dotted),  for  comparison.  The 
constant-wavelength  condition  tends  to  destabilize  peaks  but  stabilize  troughs  relative  to 
the  constant-force  condition. 

Figure  6.  Interface  profiles  dir)  for  an  initial  axisymmetric  sinusoidal  perturbation  of 
amplitude  10-3,  where  the  disturbance  is  contained  within  a  constant  radius  ( R  —  1). 
The  central  peak  grows  more  quickly  than  the  outer  ring-shaped  peak  of  the  same  initial 
amplitude. 

Figure  7.  Linearized  growth  rate  a  (A. 4c)  as  a  function  of  dimensionless  wavenumber  e: 
(a)  for  viscosity  ratio  a  =  P2/ kLi  =  10“3  and  depth  ratios  (5  =  d2/d\  —  1,  103,  and  00;  (b) 
for  deep  fluid  (/ 3  —  00)  and  a  —  10“3,  0.1,  1,  10,  103. 

Figure  8.  Linearized  growth  regimes  for  [3  3>  1  in  the  (k.  a)  parameter  plane,  showing  the 
asymptotic  growth  rates  a  in  each  regime.  (Dashed  lines  indicate  change  of  rate-controlling 
viscosity.) 

Figure  9.  Interface  profiles  d(x)  for  a  power-law  fluid,  for  an  initial  sinusoidal  perturbation 
of  amplitude  0.1,  for  the  constant-wavelength  conditions  ( L  —  1):  (a)  power-law  exponent 
m  =  3,  at  times  t  —  0,  4,  5,  and  6  (t*  «  6.08);  (b)  m  —  7  and  t  —  0,  1200.0,  and  1296.4 
(t*  ~  1296.5).  For  higher  m  the  deforming  region  is  more  compact,  and  the  profile  tends 
toward  straight  lines  between  the  deforming  regions. 

Figure  10.  Growth  of  the  disturbance  amplitude  over  time  for  a  power-law  fluid  given  an 
initial  perturbation  amplitude  of  0.1:  constant-wavelength  L  —  1  (solid),  constant  end- 
force  F  =  1  (dashed),  and  small- amplitude  approximation  (B.5)  (dotted),  (a)  m  —  3;  (b) 
m  —  7.  For  constant  force,  the  trough  (amjn)  would  pinch  off  long  before  the  peak  (ama;r) 


becomes  infinite,  but  the  constant-wavelength  condition  makes  the  peaks  less  stable  and 
troughs  more  stable.  For  higher  to,  the  growth  abruptly  becomes  rapid. 
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Figures  for  “Buoyant  Instability  of  a  Viscous  Film  Over  a  Passive  Fluid” 
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Figure  1:  Problem  description.  Between  two  horizontal,  shear-free  planes,  a  layer  of  fluid  1  of  uniform 
thickness  di  is  initially  in  unstable  equilibrium  over  a  layer  of  fluid  2  of  thickness  da-  The  upper  layer 
is  denser  (pi  >  pa),  much  thinner  (di  da))  and  much  more  viscous  (/q  pa)-  After  the  interface  is 
perturbed,  its  position  is  given  by  S(x,  y,  t). 
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Figure  2:  Balance  of  horizontal  forces,  reduced  by  the  hydrostatic  pressure  gradient  of  fluid  2,  on  a  differential 
column  of  fluid  1.  Because  there  is  no  shear  above  or  below  and  no  (reduced)  pressure  below,  the  reduced 
horizontal  force  tensor  acting  on  the  vertical  surfaces  of  the  column  has  zero  divergence. 
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Figure  3:  Interface  profiles  S(x)  at  various  times  (t  =  2.5.  3.0,  3.5.  3.7)  following  the  evolution  of  one 
wavelength  of  an  initial  sinusoidal  perturbation  of  amplitude  10~3,  for  the  constant-end-force  condition 
(F  =  1).  (Here  the  peak  reaches  infinite  thickness  at  £*  =  3.80.)  As  nonlinear  effects  become  important,  the 
peak  sharpens  while  the  troughs  flatten  and  widen;  also  the  wavelength  increases  over  time. 


Figure  4:  Interface  profiles  for  the  constant-wavelength  end  condition  (L  =  1).  Similar  to  Figure  3.  except 
that  the  troughs  become  thin  more  slowly,  the  peaks  reach  infinite  thickness  slightly  sooner  (at  t*  =  3.74). 
and  the  wavelength  remains  constant. 
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Figure  5:  Growth  of  the  disturbance  amplitude  over  time  for  an  initial  perturbation  amplitude  of  10~3. 
Both  the  peak  amplitude  am„x  =  Smax  —  1  and  the  trough  amplitude  =  1  —  Smin  are  shown:  constant- 
wavelength  condition  (solid),  constant  end-force  (dashed),  and  linearized  (exponential)  growth  (dotted),  for 
comparison.  The  constant-wavelength  condition  tends  to  destabilize  peaks  but  stabilize  troughs  relative  to 
the  constant-force  condition. 


4 


Figure  6:  Interface  profiles  S(r)  for  an  initial  axisynnnetric  sinusoidal  perturbation  of  amplitude  10_‘\  where 
the  disturbance  is  contained  within  a  constant  radius  (B  =  1).  The  central  peak  grows  more  quickly  than 
the  outer  ring-shaped  peak  of  filename  initial  amplitude. 
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Wavenumber:  e 


Figure  7:  Linearized  growth  rate  a  (A  4c)  as  a  function  of  dimensionless  wavenumber  e:  (a)  for  viscosity 
ratio  a  =  //•> / // 1  =  10~3  and  depth  ratios  3  =  d%/d\  =  1.  103,  and  oo. 
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Figure  7:  (b)  for  deep  fluid  (;3  =  oo)  and  a  =  10  3.  0.1,  1.  10.  103. 
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Figure  9:  Interface  profiles  S(x)  for  a  power-law  fluid,  for  an  initial  sinusoidal  perturbation  of  amplitude  0.1, 
for  the  constant- wavelength  conditions  (L  =  1):  (a)  power-law  exponent  m  =  3.  at  times  t  =  0.  4,  5.  and  6 
(4  «  6.08). 


Figure  9:  (b)  m  =  7  and  t  =  0.  1200.0.  and  1296.4  (£*  «  1296.5).  For  higher  m  the  deforming  region  is  more 
compact,  and  the  profile  tends  toward  straight  lines  between  the  deforming  regions. 
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Figure  10:  Growth  of  the  disturbance  amplitude  over  time  for  a  power-law  fluid  given  an  initial  perturbation 
amplitude  of  0.1:  cons f ant-' wavelength  L  =  1  (solid),  constant  end-force  F  =  1  (dashed),  and  small-amplitude 
approximation  (B  5)  (dotted).  For  Constant  force,  the  trough  (a„„:„)  would  pinch  off  long  before  the  peak 
(amax)  becomes  infinite,  but  the  constant-wavelength  condition  makes  the  peaks  less  stable  and  troughs 
more  stable.  For  higher  m.  the  growth  abruptly  becomes  rapid,  (a)  m  =  3. 
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